Exploratory analysis
Jeremy Wicquart
2020-03-05
Defining functions and packages
# 1. Source functions ----
source("functions/graphical_par.R")
source("functions/theme_graph.R")
source("functions/allometry.R")
source("functions/firstup.R")
# 2. Recquired packages ----
library(plyr)
library(tidyverse)
library(kableExtra)
library(data.table)
library(formattable)
library(rfishbase)
# Set theme_perso() as the default ggplot theme
theme_set(theme_graph())
# 3. Load data ----
data_complete <- read.csv("./../data/size_at_age_coral_reef_fishes_data.csv")Summary
## Family Genus Species
## Serranidae :1636 Epinephelus : 750 Cephalopholis argus : 518
## Acanthuridae :1225 Cephalopholis: 638 Epinephelus merra : 381
## Pomacentridae: 615 Acanthurus : 489 Myripristis berndti : 364
## Holocentridae: 522 Ctenochaetus : 429 Ctenochaetus striatus: 297
## Apogonidae : 381 Myripristis : 364 Centropyge flavissima: 251
## Balistidae : 332 Centropyge : 296 Plectropomus laevis : 248
## (Other) :1658 (Other) :3403 (Other) :4310
## ID Agei Ri Agecpt
## GAM18_B099 : 31 Min. : 0.000 Min. :0.0080 Min. : 0.000
## AC_AC_SC_03_16_001: 28 1st Qu.: 1.000 1st Qu.:0.4035 1st Qu.: 6.000
## GAM18_B007 : 28 Median : 4.000 Median :0.6062 Median : 9.000
## GAM18_B123 : 26 Mean : 4.868 Mean :0.6971 Mean : 9.738
## GAM18_A075 : 25 3rd Qu.: 7.000 3rd Qu.:0.8798 3rd Qu.:13.000
## GAM18_A168 : 25 Max. :30.000 Max. :3.7845 Max. :30.000
## (Other) :6206 NA's :377
## Rcpt Lcpt L0p R0p
## Min. :0.1517 Min. : 28.11 Min. :1.450 Min. :0.0080
## 1st Qu.:0.6445 1st Qu.:145.00 1st Qu.:1.700 1st Qu.:0.0275
## Median :0.8778 Median :206.13 Median :1.800 Median :0.0391
## Mean :1.0064 Mean :233.01 Mean :2.015 Mean :0.0436
## 3rd Qu.:1.1700 3rd Qu.:278.29 3rd Qu.:1.950 3rd Qu.:0.0544
## Max. :3.8591 Max. :984.69 Max. :4.250 Max. :0.1360
## NA's :2780
## Li_sp_m Li_sp_sd Li_sploc_m Li_sploc_sd
## Min. : 1.45 Min. : 0.0000 Min. : 1.45 Min. : 0.000
## 1st Qu.: 78.99 1st Qu.: 0.9165 1st Qu.: 74.78 1st Qu.: 1.155
## Median :146.90 Median : 2.8406 Median :145.00 Median : 3.936
## Mean :158.36 Mean : 4.3501 Mean :155.46 Mean : 5.560
## 3rd Qu.:207.97 3rd Qu.: 6.2595 3rd Qu.:205.75 3rd Qu.: 8.137
## Max. :948.14 Max. :50.2703 Max. :945.61 Max. :48.867
## NA's :398 NA's :398 NA's :745 NA's :745
## Bodymass Location
## Min. : 0.4 Gambiers :2004
## 1st Qu.: 58.0 Hao : 651
## Median : 141.3 Manuae : 113
## Mean : 530.1 Marquesas:1187
## 3rd Qu.: 313.8 Moorea :2112
## Max. :12950.0 Tuamotu : 302
## NA's :603
## Observer
## Fabien Morat :4359
## Guillemette de Synéty and Jérémy Wicquart:1368
## Jean Bienvenu : 642
##
##
##
## ## [1] 848## [1] 51## [1] 16Individual by species
data_complete %>%
group_by(Species) %>%
summarise(n = length(unique(ID)),
min_lencap = round(min(Lcpt, na.rm = TRUE), 0),
max_lencap = round(max(Lcpt, na.rm = TRUE), 0)) %>%
kable(., col.names = c("Species", "n", "Min length (TL, mm)", "Max length (TL, mm)")) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
column_spec(1, italic = T)| Species | n | Min length (TL, mm) | Max length (TL, mm) |
|---|---|---|---|
| Abudefduf sexfasciatus | 16 | 51 | 172 |
| Acanthurus achilles | 10 | 189 | 246 |
| Acanthurus lineatus | 8 | 123 | 362 |
| Acanthurus nigricans | 8 | 150 | 210 |
| Acanthurus pyroferus | 18 | 51 | 230 |
| Acanthurus triostegus | 18 | 147 | 194 |
| Balistapus undulatus | 22 | 41 | 283 |
| Caranx melampygus | 7 | 368 | 712 |
| Centropyge bispinosa | 9 | 28 | 80 |
| Centropyge flavissima | 28 | 47 | 150 |
| Cephalopholis argus | 41 | 120 | 450 |
| Cephalopholis urodeta | 8 | 150 | 214 |
| Chaetodon citrinellus | 13 | 55 | 101 |
| Chaetodon ornatissimus | 10 | 122 | 175 |
| Cheilinus chlorourus | 8 | 96 | 299 |
| Chlorurus spilurus | 34 | 164 | 344 |
| Chromis iomelas | 20 | 39 | 61 |
| Chromis viridis | 6 | 115 | 148 |
| Ctenochaetus marginatus | 15 | 99 | 270 |
| Ctenochaetus striatus | 26 | 153 | 218 |
| Dascyllus aruanus | 20 | 51 | 71 |
| Dascyllus flavicaudus | 8 | 89 | 104 |
| Epibulus insidiator | 17 | 157 | 350 |
| Epinephelus fasciatus | 10 | 123 | 256 |
| Epinephelus hexagonatus | 16 | 139 | 268 |
| Epinephelus merra | 46 | 108 | 261 |
| Epinephelus polyphekadion | 14 | 252 | 514 |
| Gnathodentex aureolineatus | 8 | 197 | 255 |
| Gymnosarda unicolor | 4 | 499 | 985 |
| Halichoeres trimaculatus | 13 | 64 | 186 |
| Lutjanus fulvus | 12 | 195 | 289 |
| Lutjanus gibbus | 17 | 113 | 155 |
| Lutjanus kasmira | 37 | 70 | 326 |
| Monotaxis grandoculis | 9 | 200 | 505 |
| Mulloidichthys flavolineatus | 14 | 257 | 311 |
| Myripristis berndti | 30 | 147 | 279 |
| Naso lituratus | 17 | 228 | 432 |
| Naso unicornis | 4 | 488 | 613 |
| Odonus niger | 18 | 176 | 350 |
| Ostorhinchus angustatus | 15 | 46 | 83 |
| Ostorhinchus apogonoides | 30 | 52 | 96 |
| Parupeneus barberinus | 16 | 82 | 496 |
| Plectropomus laevis | 31 | 356 | 957 |
| Pristiapogon taeniopterus | 32 | 51 | 104 |
| Sargocentron microstoma | 23 | 62 | 204 |
| Scarus psittacus | 8 | 195 | 311 |
| Siganus argenteus | 11 | 236 | 388 |
| Siganus spinus | 4 | 199 | 254 |
| Stegastes albifasciatus | 7 | 43 | 94 |
| Stegastes nigricans | 14 | 45 | 138 |
| Zebrasoma scopas | 18 | 110 | 166 |
Individuals by location
data_complete %>%
group_by(Species, Location) %>%
summarise(n = length(unique(ID))) %>%
pivot_wider(names_from = Location, values_from = n) %>%
mutate(Gambiers = cell_spec(Gambiers, "html", color = ifelse(is.na(Gambiers), "white", "#446CB3")),
Hao = cell_spec(Hao, "html", color = ifelse(is.na(Hao), "white", "#446CB3")),
Marquesas = cell_spec(Marquesas, "html", color = ifelse(is.na(Marquesas), "white", "#446CB3")),
Moorea = cell_spec(Moorea, "html", color = ifelse(is.na(Moorea), "white", "#446CB3")),
Manuae = cell_spec(Manuae, "html", color = ifelse(is.na(Manuae), "white", "#446CB3")),
Tuamotu = cell_spec(Tuamotu, "html", color = ifelse(is.na(Tuamotu), "white", "#446CB3"))) %>%
kable(., format = "html", escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
column_spec(1, italic = T)| Species | Gambiers | Moorea | Manuae | Marquesas | Tuamotu | Hao |
|---|---|---|---|---|---|---|
| Abudefduf sexfasciatus | 11 | 5 | NA | NA | NA | NA |
| Acanthurus achilles | NA | NA | 10 | NA | NA | NA |
| Acanthurus lineatus | NA | NA | NA | 8 | NA | NA |
| Acanthurus nigricans | NA | NA | NA | 8 | NA | NA |
| Acanthurus pyroferus | NA | NA | NA | 18 | NA | NA |
| Acanthurus triostegus | 3 | 11 | NA | 4 | NA | NA |
| Balistapus undulatus | NA | 14 | NA | 8 | NA | NA |
| Caranx melampygus | NA | 5 | NA | NA | 2 | NA |
| Centropyge bispinosa | NA | 9 | NA | NA | NA | NA |
| Centropyge flavissima | 9 | 7 | NA | 12 | NA | NA |
| Cephalopholis argus | 6 | 10 | NA | NA | 10 | 15 |
| Cephalopholis urodeta | 8 | NA | NA | NA | NA | NA |
| Chaetodon citrinellus | NA | 13 | NA | NA | NA | NA |
| Chaetodon ornatissimus | 4 | NA | NA | 6 | NA | NA |
| Cheilinus chlorourus | 2 | 6 | NA | NA | NA | NA |
| Chlorurus spilurus | 13 | 19 | NA | NA | 2 | NA |
| Chromis iomelas | NA | 20 | NA | NA | NA | NA |
| Chromis viridis | 6 | NA | NA | NA | NA | NA |
| Ctenochaetus marginatus | NA | NA | NA | 15 | NA | NA |
| Ctenochaetus striatus | 11 | 15 | NA | NA | NA | NA |
| Dascyllus aruanus | NA | 20 | NA | NA | NA | NA |
| Dascyllus flavicaudus | 8 | NA | NA | NA | NA | NA |
| Epibulus insidiator | 12 | 5 | NA | NA | NA | NA |
| Epinephelus fasciatus | NA | NA | NA | 10 | NA | NA |
| Epinephelus hexagonatus | NA | NA | NA | NA | NA | 16 |
| Epinephelus merra | 13 | 20 | NA | NA | NA | 13 |
| Epinephelus polyphekadion | NA | NA | NA | NA | NA | 14 |
| Gnathodentex aureolineatus | 8 | NA | NA | NA | NA | NA |
| Gymnosarda unicolor | 4 | NA | NA | NA | NA | NA |
| Halichoeres trimaculatus | 7 | 6 | NA | NA | NA | NA |
| Lutjanus fulvus | NA | 6 | NA | NA | 6 | NA |
| Lutjanus gibbus | NA | NA | NA | 17 | NA | NA |
| Lutjanus kasmira | 2 | NA | NA | 35 | NA | NA |
| Monotaxis grandoculis | 5 | 3 | NA | NA | 1 | NA |
| Mulloidichthys flavolineatus | NA | 14 | NA | NA | NA | NA |
| Myripristis berndti | 7 | 15 | NA | 8 | NA | NA |
| Naso lituratus | 11 | 6 | NA | NA | NA | NA |
| Naso unicornis | 4 | NA | NA | NA | NA | NA |
| Odonus niger | NA | 13 | NA | 5 | NA | NA |
| Ostorhinchus angustatus | NA | 15 | NA | NA | NA | NA |
| Ostorhinchus apogonoides | NA | NA | NA | 30 | NA | NA |
| Parupeneus barberinus | 5 | 4 | NA | 7 | NA | NA |
| Plectropomus laevis | 12 | NA | NA | NA | NA | 19 |
| Pristiapogon taeniopterus | NA | NA | NA | 32 | NA | NA |
| Sargocentron microstoma | NA | 17 | NA | 6 | NA | NA |
| Scarus psittacus | NA | 3 | NA | NA | 5 | NA |
| Siganus argenteus | 6 | 2 | NA | NA | 3 | NA |
| Siganus spinus | NA | 4 | NA | NA | NA | NA |
| Stegastes albifasciatus | NA | 7 | NA | NA | NA | NA |
| Stegastes nigricans | 3 | 11 | NA | NA | NA | NA |
| Zebrasoma scopas | 5 | 13 | NA | NA | NA | NA |
Individuals with back-calculation
# 1. Create a dataset without NA ----
data_without_na <- data_complete %>%
filter(!is.na(Li_sp_m))
# 2. Get the number of species with back-calculated size-at-age ----
length(unique(data_without_na$Species))## [1] 45# 3. Get the number of individuals with back-calculated size-at-age ----
length(unique(data_without_na$ID))## [1] 709# 4. Get the number of individuals with back-calculated size-at-age ----
length(unique(data_without_na$Family))## [1] 16Agei vs Ri
ggplot(data_complete, aes(x = Agei, y = Ri)) +
geom_point(color = col_color_graph, fill = col_fill_graph, size = 1, shape = 21)+
theme(strip.text.x = element_text(size = 8, face = "italic"),
strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
plot.title = element_text(colour = col_color_graph),
plot.subtitle = element_text(colour = "black")) +
facet_wrap(~Species, scales = "free", ncol = 5) +
labs(x = "Age (years)", y = "Radius at age i (mm)", title = "Check for outliers", subtitle = "agei vs radi")
Agei vs Li (species)
ggplot(data_complete, aes(x = Agei, y = Li_sp_m)) +
geom_point(color = col_color_graph, fill = col_fill_graph, size = 1, shape = 21)+
theme(strip.text.x = element_text(size = 8, face = "italic"),
strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
plot.title = element_text(colour = col_color_graph),
plot.subtitle = element_text(colour = "black"),
strip.text.y = element_text(angle = 360)) +
facet_grid(Species~Location, scales = "free") +
labs(x = "Age (years)", y = "Length at age i (mm)", title = "Check for outliers", subtitle = "agei vs li")
Agei vs Li (species and location)
ggplot(data_complete, aes(x = Agei, y = Li_sploc_m)) +
geom_point(color = col_color_graph, fill = col_fill_graph, size = 1, shape = 21)+
theme(strip.text.x = element_text(size = 8, face = "italic"),
strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
plot.title = element_text(colour = col_color_graph),
plot.subtitle = element_text(colour = "black"),
strip.text.y = element_text(angle = 360)) +
facet_grid(Species~Location, scales = "free") +
labs(x = "Age (years)", y = "Length at age i (mm)", title = "Check for outliers", subtitle = "agei vs li")
Agei vs Li (comparison)
data_complete %>%
select(-Li_sp_sd, -Li_sploc_sd) %>%
pivot_longer("Li_sp_m":"Li_sploc_m", names_to = "Method", values_to = "Li") %>%
mutate(Method = str_replace_all(Method, c("Li_sp_m" = "Species",
"Li_sploc_m" = "Species x Location"))) %>%
ggplot(., aes(x = Agei, y = Li, fill = Method)) +
geom_point(color = col_color_graph, size = 1, shape = 21, alpha = 0.5)+
theme(strip.text.x = element_text(size = 8, face = "italic"),
strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
plot.title = element_text(colour = col_color_graph),
plot.subtitle = element_text(colour = "black"),
strip.text.y = element_text(angle = 360),
legend.position = "top") +
facet_grid(Species~Location, scales = "free") +
labs(x = "Age (years)", y = "Length at age i (mm)")
Bodymass vs Lcpt
# 1. Remove NA from data ----
data_complete_na_rm <- data_complete %>%
filter(!is.na(Bodymass) & !is.na(Lcpt))
# 2. Get the parameters of the allometric regression for each species ----
coef_allo <- plyr::ddply(data_complete_na_rm, .(Species), function(x) allometry(data = x, xreg = "Lcpt", yreg = "Bodymass", graph = FALSE))
# 3. Generate predicted values ----
pred_values <- plyr::ddply(data_complete_na_rm, .(Species), function(x) seq(0, max(x$Lcpt), length.out = 100)) %>% # Create new x
gather(., pos, xnew, -1) %>% # Gather multiple columns into two columns
select(Species, xnew) %>% # Delete useless column
arrange(Species) %>% # Sort DF by species
merge(., coef_allo, by = "Species") %>% # Merge the DF with parameters of allometric regression
mutate(ynew = a*xnew^b) # Calculate new y
# 4. Plot ----
ggplot() +
geom_line(data = pred_values, aes(x = xnew, y = ynew), color = col_color_graph) +
geom_point(data = data_complete, aes(x = Lcpt, y = Bodymass), color = col_color_graph, fill = col_fill_graph, size = 1, shape = 21)+
theme(strip.text.x = element_text(size = 8, face = "italic"),
strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
plot.title = element_text(colour = col_color_graph),
plot.subtitle = element_text(colour = "black")) +
facet_wrap(~Species, scales = "free", ncol = 5) +
labs(x = "Length at capture (mm)", y = "Bodymass at capture (g)", title = "Check for outliers", subtitle = "lencap vs Bodymass")
Export table for Morat et al
# 1. Get the type of each variable ----
table_type <- data_complete %>%
summarise_all(class) %>%
gather(.) %>%
rename(Variable = key, Type = value) %>%
mutate(Type = firstup(Type), # First letter in uppercase
Column = 1:nrow(.)) # Add the numbering of the variable
# 2. Get the range of each numeric variable ----
table_range <- data_complete %>%
summarise_if(is.numeric, ~paste(round(min(.x, na.rm = TRUE), digits = 3),
round(max(.x, na.rm = TRUE), digits = 3),
sep = " - ")) %>%
gather(.) %>%
rename(Variable = key, Range = value)
# 4. Merge the three datasets and rename variables ----
merge(table_type, table_range, by = "Variable", all.x = TRUE) %>%
mutate(Description = NA, Unit = NA) %>%
select(Column, Variable, Description, Unit, Type, Range) %>%
arrange(Column) %>%
write.csv2(., "./../data/variables_complete_for_Morat_et_al.csv", row.names = FALSE)Number of NA by variable
data_complete %>%
summarise_all(~(sum(is.na(.)))) %>%
t(.) %>%
as.data.frame() %>%
rename("NA" = "V1")## NA
## Family 0
## Genus 0
## Species 0
## ID 0
## Agei 0
## Ri 377
## Agecpt 0
## Rcpt 0
## Lcpt 0
## L0p 0
## R0p 2780
## Li_sp_m 398
## Li_sp_sd 398
## Li_sploc_m 745
## Li_sploc_sd 745
## Bodymass 603
## Location 0
## Observer 0Difference between sp and site x sp
data_lm <- lm(Li_sploc_m ~ Li_sp_m, data = data_complete)
data_complete %>%
ggplot(., aes(x = Li_sp_m, y = Li_sploc_m)) +
geom_point() +
geom_abline(slope = data_lm$coefficients[2], intercept = data_lm$coefficients[1]) +
annotate("text", x = Inf, y = Inf,
label = paste("r² = ", round(summary(data_lm)$r.squared, 3), sep = ""),
vjust = 1, hjust = 1) +
labs(x = "Length at age i (sp)", y = "Length at age i (sp x site)")
Check lencap vs fishbase
# 1. Download fishbase max size for each species ----
fishbase_maxlength <- length_length(unique(data_complete$Species), fields = c("Species", "LengthMax", "Type")) %>%
filter(!is.na(LengthMax), Type == "TL") %>%
select(-Type) %>%
group_by(Species) %>%
filter(LengthMax == max(LengthMax, na.rm = TRUE)) %>%
mutate(LengthMax = LengthMax*10) %>%
ungroup()
# 2. Graphic observations before removing outliers ----
data_complete %>%
select(Species, Lcpt) %>%
ggplot(data = .) +
geom_point(aes(x = Species, y = Lcpt), alpha = 0.6, shape = 16, fill = "#446CB3", color = "#446CB3") +
geom_boxplot(aes(x = Species, y = Lcpt), fill = "#ECECEC", alpha = 0.5) +
geom_point(data = fishbase_maxlength, aes(x = Species, y = LengthMax), size = 1.2, color = "red", fill = "red", shape = 23) +
coord_flip() +
labs(x = NULL, y = "Total length (mm)")
Reproducibility
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
## [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
## [5] LC_TIME=French_France.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rfishbase_3.0.4 formattable_0.2.0.1 data.table_1.12.8
## [4] kableExtra_1.1.0 forcats_0.5.0 stringr_1.4.0
## [7] dplyr_0.8.4 purrr_0.3.3 readr_1.3.1
## [10] tidyr_1.0.2 tibble_2.1.3 ggplot2_3.2.1
## [13] tidyverse_1.3.0 plyr_1.8.6 RColorBrewer_1.1-2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 lubridate_1.7.4 lattice_0.20-38 assertthat_0.2.1
## [5] digest_0.6.25 R6_2.4.1 cellranger_1.1.0 backports_1.1.5
## [9] reprex_0.3.0 evaluate_0.14 httr_1.4.1 highr_0.8
## [13] pillar_1.4.3 rlang_0.4.5 curl_4.3 lazyeval_0.2.2
## [17] readxl_1.3.1 rstudioapi_0.11 rmarkdown_2.1 labeling_0.3
## [21] webshot_0.5.2 htmlwidgets_1.5.1 munsell_0.5.0 broom_0.5.5
## [25] compiler_3.6.3 modelr_0.1.6 xfun_0.12 pkgconfig_2.0.3
## [29] htmltools_0.4.0 tidyselect_1.0.0 fansi_0.4.1 viridisLite_0.3.0
## [33] crayon_1.3.4 dbplyr_1.4.2 withr_2.1.2 grid_3.6.3
## [37] nlme_3.1-144 jsonlite_1.6.1 gtable_0.3.0 lifecycle_0.1.0
## [41] DBI_1.1.0 magrittr_1.5 scales_1.1.0 cli_2.0.2
## [45] stringi_1.4.6 reshape2_1.4.3 farver_2.0.3 fs_1.3.1
## [49] xml2_1.2.2 ellipsis_0.3.0 generics_0.0.2 vctrs_0.2.3
## [53] gh_1.1.0 tools_3.6.3 glue_1.3.1 hms_0.5.3
## [57] yaml_2.2.1 colorspace_1.4-1 rvest_0.3.5 memoise_1.1.0
## [61] knitr_1.28 haven_2.2.0Jeremy WICQUART | jeremywicquart@gmail.com | 2020-03-05 17:13:34